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METHOD AND DEVICE FOR THE FOURTH-ORDER, BLIND 
IDENTIFICATION OF AN UNDER-DETERMINED MIXTURE OF SOURCES 

BACKGROUND OF THE INVENTION 

1. Field of the Invention 

The invention relates especially to a method for the learned or 
blind identification of a number of sources P that is potentially greater than 
5 or equal to the number N of sensors of the reception antenna. 

It can be used for example in the context of narrow-band multiple 
transmission. 

It is used for example in a communications network. 
It can be applied especially in the field of radio communications, 
10 space telecommunications or passive listening to these links in frequencies 
ranging for example from VLF to EHF. 

Figure 1 is a schematic drawing exemplifying an array of several 
reception sensors or receivers, each sensor receiving signals from one or 
more radio communications transmitters from different directions of arrival 
15 Each sensor receives signals from a source with a phase and 

amplitude that are dependent on the angle of incidence of the source and 
the position of the sensor. Figure 2 is a drawing exemplifying the 
parametrization of the direction of a source. This direction is parametrized 
by two angles corresponding to the azimuth angle 8 and the elevation angle 
20 A. 

2. Description of the Prior Art 

The past 15 years or so have seen the development of many 
techniques for the blind identification of signatures or source direction 
vectors, assumed to be statistically independent. These techniques have 

25 been developed in assuming a number of sources P smaller than or equal to 
the number of sensors N. These techniques have been described in the 
references [1][3][7] cited at the end of the description. However, for many 
practical applications such as HF radio communications, the number of 
sources from which signals are received by the sensors is increasing 

30 especially with the bandwidth of the receivers, and the number of sources P 
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can therefore be greater than the number of sensors N. The mixtures 
associated with the sources are then said to be under-determined. 

A certain number of methods for the blind identification of under- 
determined mixtures of narrow-band sources for networks have been 
5 developed very recently and are described in the references [2] [7-8] and 
[10]. The methods proposed in the references [2] and [7-8] make use of 
information contained in the fourth-order (FO) statistics of the signals 
received at the sensors while the method proposed in the reference [10] 
make use of the information contained in one of the characteristic functions 

10 of the signals received. However, these methods have severe limitations in 
terms of the prospects of their operational implementation. Indeed, the 
method described in the reference [2] is very difficult to implement and does 
not provide for the identification of the sources having the same kurtosis 
values (standardized fourth-order cumulant). The methods described in the 

15 references [7-8] assume that the sources are non-circular. These methods 
give unreliable results in practice. Finally, the reference method [10] has 
been developed solely for mixtures of sources with real (non-complex) 
values. 

The object of the present invention relates especially to a new 
20 method for the blind identification of an under-determined mixture of 
narrow- band sources for communications networks. The method can be 
used especially to identify up to N 2 - N + 1 sources from N identical sensors 
and up to N 2 sources with N different sensors, in assuming only that the 
sources have different tri-spectra and non-zero, same-sign kurtosis values 
25 (this hypothesis is practically always verified in the context of radio 
communications) . 

SUMMARY OF THE INVENTION 

The invention relates to a method for the fourth-order, blind 
identification of at least two sources in a system comprising a number of 
30 sources P and a number N of reception sensors receiving the observations, 
said sources having different tri-spectra, wherein the method comprises at 
least the following steps: 
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o a step for the fourth-order whitening of the observations received on 
the reception sensors in order to orthonormalize the direction vectors 
of the sources in the matrices of quadricovariance of the observations 
used, 

5 o a step for the joint diagonalizing of several whitened matrices of 
quadricovariance in order to identify the spatial signatures of the 
sources. 

The number of sources P is for example greater than the number 
of sensors N. 

10 The method can be used in a communication network. 

The method according to the invention has especially the 
following advantages: 

o it enables the identification of a number of sources greater than the 
number of sensors: 

15 o for identical sensors: N 2 - N + 1 sources having different tri- 

spectra and non-zero and same-sign kurtosis values, 
o for different sensors (arrays with polarization diversity and/or 
pattern diversity and/or coupling etc) N 2 sources having 
different tri-spectra and non-zero, same-sign kurtosis values, 
20 o the method is robust with respect to Gaussian noise, even spatially 
correlated Gaussian noise, 
o it enables the goniometry of each source identified, using a wavefront 
model attached to the signature, with a resolution potentially higher 
than that of existing methods, 
25 o it enables the identification of I (N 2 - N + 1) cyclostationary sources if 
the sensors are identical and I x N 2 cyclostationary sources if the 
sensors are different: with polarization diversity and/or pattern 
diversity and/or coupling, where I is the number of cyclical 
frequencies processed, 
30 o using a performance criterion, it enables the quantitative evaluation 
of the quality of estimation of the direction vector of each source and 
a quantitative comparison of two methods for the identification of a 
given source, 
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o using a step for the selection of the cyclical frequencies, it enables 
the processing of a number of sources greater than the number of 
sources processed by the basic method. 

BRIEF DESCRIPTION OF THE DRAWINGS 

Other features and advantages of the invention shall appear more 
clearly from the following description along with the appended figures, of 

which: 

o Figure 1 shows an example of a communication network, 
o Figure 2 shows parameters of a source, 

o Figure 3 is a functional diagram of the method according to the 
invention, 

o Figure 4 shows an example of spatial filtering, 

o Figures 5 and 6 show examples of variations of the performance 
criterion as a function of the number of samples observed, comparing 
the performance of the method with two prior art methods. 

o Figures 7 to 9 show three alternative embodiments of the method 
described in figure 3 implementing a selection of the cyclical 
frequencies. 

MORE DETAILED DESCRIPTION 

For a clear understanding of the object of the invention, the 
following example is given by way of an illustration that in no way restricts 
the scope of the invention for a radio communications network in a 
multiple-transmission, narrow-band context, with sources having different 
tri-spectra (of cumulants). 

Each sensor of the network, formed by N receivers, receives a 
mixture of P narrow-band (NB) sources which are assumed to be 
statistically independent. On this assumption, the vector of the complex 
envelopes of the signals at output of the sensors is written as follows: 

p 

xft/=X s P (t)a P + b(t) =As(t) + b(t) (1) 

P =\ 

where s p (t) is the signal of the p th source as well as the component of the 
vector s(t ), b{t ) is the assumed Gaussian noise vector with any covariance, 
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a p is the signature or the direction vector of the source and A (N x P) is 
the matrix of the vectors a p (direction vectors of the sources). 

It is an object of the invention especially to identify the direction 
vectors a p of each of the sources when, especially, the number of sources P 
5 is potentially greater than the number of sensors N. 

From this identification, it is then possible to apply techniques for 
the extraction of the sources by the spatial filtering of the observations. The 
blind extraction is aimed especially at restoring the information signals 
conveyed by the sources in not exploiting any a priori information (during 
10 normal operation) on these sources. 

Fourth-order statistics 

The method according to the invention makes use of the fourth-order 
statistics of the observations corresponding to the time-domain averages 
Qa(ti,T2,T3) =< Qx(n,T2,T3)(t)>, on an infinite horizon of observation, of certain 
15 matrices of quadricovariance, £MTi,T 2 ,T 3 )(f), sized (N^xN 2 ). The elements, 
Qx(ti,T2,T3)[i, j, K l\(t) of these matrices are, for example, defined by the 
relationship: 

£?.v(t,,t 2 ,t 3 )[z, j, kj][t] = Cum(xj(t),Xjit-T { )\xk(t'T 2 )\ x/(/ -t 3 )) (2) 
20 where * is the conjugate complex symbol, x((t) is the I th component of the 
vector x(t) , <.> is the operation of time-domain averaging on an infinite 
horizon of observation and (ti,T2,t 3 ) is a triplet of delays. Assuming that 
CMti,T2,T3)[z, j, k, /] is the element [N(i -1) + j, N(k - 1) + Z] of the matrix 
CMn,T2,T3), assuming that the noise is Gaussian and using the expression 
25 (1) in the expression (2), the matrix of quadricovariance CMti,t 2 ,t 3 ) is written 
as follows: 

g,( r It t 2t t 3 ) = (A ®A*) Q s ( r h t 2 , r>) (A ®A*f (3) 

where Qs(ti,T2,T3) is the averaged matrix of quadricovariance of s{t) with a 
dimension (i^xP 2 ), A=[ a\ ... a P ], ® is the Kronecker product and H 
30 designates the transpose and conjugate. Assuming statistically independent 
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sources, the matrix Qs(xi,X2,x 3 ) is formed by at least P 4 - P zeros and the 

expression (3) is simplified as follows: 

p 

£?.r(X|,T 2 ,x 3 ) = J] c / Xt 1 ,t 2 ,t 3 ) {a p ® a/) (a p ® a/f (4a) 

p=\ 

= Aq Q(X„X 2 ,X 3 ) V (4b) 

5 where Ag is a matrix sized {ifx F) defined by Aq = [(ai®ai*), (a p ®a p *)], 

C s (ti,i2,i3) is a diagonal matrix sized (P x F) defined by C s (xi,X2,X3) = 

diag[ci(Ti,T2,T3),..., Cp(n,T2,T3)] and where c p (ti, 12,13) is defined by: 

c^(t,,X2,t 3 ) = <Cum(s p (tl Sp(/-X,)\ s p (t-z 2 )\ 5p(M 3 ))> (5) 

The expression (4b) has an algebraic structure similar to that of 

10 the correlation matrix of the observations used in the algorithm SOBI 

described in the reference [1]. The notation used here below will be Q x = 

Q x (0, 0, 0), c p = c p (0 } 0, 0), C s = C s (0, 0, 0) in order to deduce the relationship 

(4b) therefrom: 

Q x = Aq C s AqH (6) 

15 It is assumed here below that the number of sources P is such 

that P < rf y that the matrix Aq is of full rank, that the averaged cumulants 
c p , 1 < p < P, are non-zero (non-Gaussian sources) and same-sign 
cumulants and that, for any pair (i, j) of sources, there is at least one triplet 
of delays (11,12,13) such that | xi | + | t 2 | + 1 13 1 * 0 and 

20 C I (Xl,X2,X 3 ) / \d\ * C,(X1,X2,X 3 ) / |Q| (7) 

Fourth-order whitening step 

The first step of the method according to the invention, called 
FOBIUM, consists of the orthonormalization, in the matrix of 
25 quadricovariance Q x of the expression (6), of the columns of the matrix A Qy 
considered to be virtual direction vectors of the sources for the array of 
sensors considered. To this end, the method considers the eigen-element 
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decomposition of the P rank hermitian matrix Q x given by 

Q x = E X A X E X H (8) 

where A* is the real longitudinal diagonal, with a dimension (P x P), of the P 

non-zero eigenvalues of Q x , and E x is the matrix sized (N*xF) of the 

5 associated, orthonormal eigenvectors. For a full-rank matrix Aq, it can be 

shown that there is equivalence between assuming that the kurtosis values 

of the sources have a same sign e (e = ±1) and assuming that the 

eigenvalues of A x also have a same sign e. In this context, it is possible to 

build the following whitening matrix T sized (Px iV 2 ): 
10 T = (A x )~ m E X H (9) 

The whitening matrix sized (P x pf) is defined from the square 
root of the real diagonal matrix sized (P x P) of the P non-zero eigenvalues 
of the matrix of quadricovariance and of the transpose of the matrix of the 
associated eigenvectors with a dimension (Px N 2 ] 

15 where (A*) 1/2 is the inverse of the square root of A x . From the expressions (6) 

and (8), it is deduced that: 

sTQ x f = TA Q (sC s )A Q H T H = I, (10) 

where I P is the identity matrix with a dimension (P x P) and where eC s = 

diag[|ci|,..., \Cp\]. This expression shows that the matrix TA Q (sC s ) 1/2 with a 

20 dimension (P x P\ is a unitary matrix U. It is deduced from this that: 

TA Q = U{sCs) m (11) 

Fourth-order identification step 

From the expressions (4b) and (11), it is deduced that: 

25 ^(T.,I2,X3) = ^,(T I ,T 2 ,T 3 )7^=C/(^Q)" ,/2 C J (X I ,T 23 T3) (sCs)'^ if (12) 

where W(x\ y 12,13) is the matrix of quadricovariance whitened at the fourth 
order by the matrix Q x . This expression shows that the unitary matrix U 
diagonalizes the matrices T Qx(ti,t 2 ,t 3 ) and that the associated 
eigenvalues are the diagonal terms of the diagonal matrix (eC s ) m C s (ti,T2,t 3 ) 
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(eC s ) . For a given triplet of delays (11,12,13), the matrix U is unique, give or 

take one permutation and one unit diagonal matrix, when the elements of 

the matrix (eC s ) 1 ~ C s (ti,T2,T3) (eC s ) ^all are different. If not, the method uses 

a set of K triplets (T\ k y X2 k ,t$ k ) , l< k < K y defined as follows: for all pairs of 

5 sources (i, j), there is at least one triplet (xi k f X2 k f X3 k ) y such that the condition 

of the equation (7) is verified. In these conditions, the unitary matrix U is 

the only matrix U so i which, to the nearest permutation and to the nearest 

unit diagonal matrix, jointly diagonalizes the K matrices TxQ x (x\ k ) X2 k ,X3 k )xT yi . 

Consequently, the matrix U S ot, which resolves the above problem, is written 

10 as a function of the unitary matrix U as follows: 

U sol =U An (13) 

where A and n are respectively the unit diagonal matrix and the 

permutation matrix referred to here above. From the equations (11) and 

(13), it is possible to deduce the matrix A Q to the nearest unit diagonal 

15 matrix and to the nearest permutation, which is expressed by: 

fu so] =[b^..b P )= E x V^soi =A Q (£C s ) a ATl (14) 

where f is the pseudo-inverse of the matrix T. Each column, bi (1< I < F) } of 
the matrix f U so i corresponds to one of the vectors |i q \c q \ ]/2 (a q ®a q *} ) 1< q < 
P, where m is a complex scalar value such that |^ q | = 1. Consequently, in 
20 converting each column b\ of the matrix f U so t into a matrix B\ with a 
dimension (N x TV) such that j\ = bi((i -1)N + j) (1 < i, j < N), it is deduced 
therefrom that: 

1/2 u 

B f =[i (l \c t/ \ a q a q H pour (1< (15) 
The matrix B ( is built from the vector b[ and depends on a 
25 complex scalar value, the square root of the cumulant and the direction 
vector of the q th source and of its conjugate. 

In this context, the direction vector Oq of the cp* source is associated with 
the eigenvector of B\ associated with the highest eigenvalue. 
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Summary of the principle of the invention 

In brief, the different steps of the method according to the 
invention include at least the following steps: 

for L vector observations received in the course of the time: x(lT e ) (1< I < L), 
5 where T e is the sampling period. 

Estimation 

Step 1: The estimation, through Q A X) of the matrix of quadricovariance Q x , 
from the L observations x(lT e ), using a non-skewed and asymptotically 
consistent estimator. Depending on the nature of the sources, the estimator 
10 is adapted as follows: 

• Stationary and centered case: empirical estimator used in the 
reference [3]. 

• Cyclostationary and centered case: estimator implemented in the 
reference [10]. 

15 • Cyclostationary and non-centered case: estimator implemented in the 
reference [11]. 
Whitening 

Step 2: The eigen-element decomposition of the estimated matrix of 
quadricovariance Q A Xl estimating the number of sources P and restricting 
20 this eigenvalue decomposition to the P main components: Q; A X « E A X A; A X 
E] A X H , where A; A * is the diagonal matrix containing the P highest modulus 
eigenvalues and E A X is the matrix containing the associated eigenvectors. 
Step 3: The building of the whitening matrix: T A = (A' A X )' U2 E A X ". 

25 Selection of the triplets 

Step 4: The selection of K triplets of delays (ti*,T2 /c ,t 3 ' c ) where 

|T,M + |T2M + |T3M*0. 

Estimation 
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Step 5: The estimation, through Q; A J c(Ti fc ,T2 fc ,T 3 * : ), of K matrices of 
quadricovariance Q x (Ti fc ,T2 fc ,T3 fc ). As in the step 1, this estimation depends 
especially on the assumptions made on the observations: 

o Stationary and centered case: empirical estimator used in the 
5 reference [3]. 

o Cyclostationary and centered case: estimator implemented in the 
reference [10]. 

o Cyclostationary and non-centered case: estimator implemented in the 
reference [11]. 
10 Identification 

Step 6: The computation of the matrices T A Q A x (xi k ,X2 k ,X3 k ) T A " and the 
estimation, by U; A S ot, of the unitary matrix U so i through the joint 
diagonalizing of the K matrices T A Q A x (x\ k y X2 k 7 T3 k ) T Ah . 

Step 7: The computation of T; A V ; A S0 /=[b; A i...b ;%] and the building of the 

15 matrices B; A / sized (TV x TV). 

Step 8: The estimation, by a; A P , of the signatures Oq (1< q < P) of the P 
sources in carrying out a decomposition into elements on each matrix B; A L 
Applications 

At the end of the step 8, the method has identified the direction 
20 vectors of P non-Gaussian sources having different tri-spectra with same- 
sign kurtosis values. P<N 2 and P may reach N 2 - N + 1 or N 2 depending on 
the type of sensors used. 

Using this information, the method may implement a method of 
goniometry or a spatial filtering of antennas. 
25 A method of goniometry can be used to determine the direction of 

arrival of the sources and more precisely the azimuth angle 0 m for ID 
goniometry and azimuth and elevation angles (9 m , A m ) for 2D goniometry. 

Figure 4 represents a spatial filtering of antennas for spatial 
filtering structures. It enables especially the optimizing of reception from 



one or all the sources present by the spatial filtering of the observations. 
When several sources are of interest to the receiver, we speak of source 
separation techniques. When no a priori information on the sources is 
exploited, we speak of blind techniques. 

Verification of the quality of the estimates 

According to one alternative embodiment, the method comprises 
a step of qualitative evaluation, for each source, of the quality of 
identification of the associated direction vector. 

This new criterion enables the intrinsic comparison of two 
methods of identification for the restitution of the signature of a particular 
source. This criterion, for the identification problem, is an extension of the 
one proposed in [5] for extraction. It is defined by the P-uplet : 

D(Aj) = (a,,a 2 , , a P ) (16) 

where 

a^min^a,)] (17) 

and where d(u,v) is the pseudo-distance between the vectors u and v, such 
that: 

|„"vf 

v) = ' - Rqk) (,8) 

In the simulations of figures 5 and 6, there are P=6 statistically 
independent sources received on a circular array of N=3 sensors having a 
radius r such that r/A.=0.55 (X: wavelength). The six sources are non-filtered 
QPSK sources having a signal-to-noise ratio of 20dB with a symbol period T 
= 4T e , where T e is the sampling period. 

The incidence values of the sources are such that 0i=2.16°, 
&=25.2°, ft=50°, 04=272.16°, ft=315.36°, 0 6 =336.96° and the associated 
carrier frequencies verify A/\ T e =0, Af 2 7W/2, A/ 3 7>l/3, A/ 4 7>l/5, A/ 5 
T e =l/7 and Af 6 7W/11. The JADE [3], SOBI [1] and FOBIUM method 
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according to the invention are applied and the performance values a q for 
q=1...6 are evaluated after an averaging operation on 1000 results. For the 
FOBIUM method, we choose K=4 triplets of delays (n k i T2 k ,T3 k ) where r\ k =kT e 
and T2 k =Tz k =0. 

In the above assumptions, figure 5 shows the variation in <x 2 
(performance of the second source) resulting from the JADE (b), SOBI (c) 
and FOBIUM (a) separators as a function of the number L of samples. The 
curves show firstly that the JADE and SOBI methods present difficulties in 
identifying the direction vector of the second source in an under-determined 
mixture context and that, secondly, that the FOBIUM method performs very 
well. 

Figure 6 gives a view, in the same context, of the variations of all 
the a p (1 < p < 6) values resulting from the FOBIUM method as a function of 
L. The curve (index p) is associated with the p th source. It can be seen that 
all the coefficients a p converge towards zero and that, asymptomatically, the 
direction vectors are perfectly identified. 
Variants of the cyclical FOBIUM method 

Figures 7 and 8 show two examples of variants of the method 
according to the invention known as the cyclical FOBIUM method. 

The idea lies especially in introducing selectivity by the cyclical 
frequencies into the method presented here above and is aimed especially at 
the blind identification, with greater processing capacity, of under- 
determined mixtures of cyclo stationary sources. 

The major difference between the steps 1 to 8 explained here 
above and this variant is the implementation of a step for the cyclical 
isolation of the sources by fourth-order discrimination according to their 
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cyclical frequencies. It is thus possible to separately identify the sources 
associated with a same fourth-order cyclical parameter without being 
disturbed by the other sources processed separately. 

The two variants shown in figures 7 and 8 can be implemented 
5 by reiterating the process of cyclical isolation on the "other sources" with 
other cyclical parameters. The process of cyclical isolation can be applied 
several times in a third version illustrated in figure 9. 

Fourth-order cyclical statistics 

10 The fourth-order cyclical statistics of the observations or sensor 

signals used are characterized by the matrices of cyclical qu ad rico variance 
£>/(a,xi,T2,T3) with a dimension (N 2 x N 2 ) where the elements Q/(a,xi,X2,X3)[z, j, 
/c, /], are defined by: 

gvWi^Ts) k, I] = <Cum(x/(0, xj(t -t,)\ xj&t -x 2 )\ x\(t -t 3 )) exp(-j27ia/)> 
15 0, 2 (a,T,,T 2 ,T 3 ) k, I] = <Cum(x/(/), xj(t -t,) * , x&t -t 2 ) , x[(t -t 3 )) exp(-j2ra0> 

a 3 (a,T,,T2,T3) [ij, K i] = <Cum(x/(/), xj(t -x,) , x&t -x 2 ) , x{(t -x 3 )) exp(-j2rca/)> (19) 

It can be seen that g/(a,x h x 2 ,x 3 ) is associated with the e th fourth-order 
moment. In stating that £?/(a,X|,x 2 ,x 3 )[/,y, k t I] is the element [N(i -1) + j, N(k- 1) 
+ /] of the matrix (?/(a,X|,x 2 ,x 3 ) and assuming that the noise is Gaussian, the 
20 matrix £?/(a,X|,x 2 ,x 3 ), is written as follows, in using (1) and (19): 
&W/.r ;> f>) = (A®A*)Q s \a y T h T 2 ,T s )(A®A*)" 
a 2 (a, v h rji rj) = (A ®A* ) Q s \a, r h r 2 , z 3 ) (A ®A) J 
g/(a, r ; , t : , t s ) = (A®A) Q s \a, r,, r 2 , z>) (A ®A) T (20) 

Where Q s c (a,xi,x 2 ,x 3 ) is a cyclical matrix of quadricovariance s(t) with a 
25 dimension (P 2 xP 2 ), ® is the Kronecker product and T designates the 
transpose. On the assumption of statistically independent sources, the 
matrix Qs*(a,xi,x 2 ,x 3 ) is formed by at least P 4 - P zeros and the expression 
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(20) is simplified as follows: 

p 

£?,Wi,t 2 ,t 3 ) = £ c /'( a ^i^3)(^fl/)K®«/) H = A Q C s \a,x { ,x 2i x 3 ) A q h 

p=\ 

p 

£?;W,,t 2 ,T3)= £ c^a^i^.^OCa^^^Cfl/,®^,) 7 = C s 2 (a,x u x 2 ,xi) 5<? T 

e,Wi,X2,i3)= X ^ 3 (a ) x l ,x 2 ,x 3 )K®^)K®fl / ,) T = fl G C s 3 (a,Tb* 2 ,T 3 ) V ( 21 > 
/>=l 

5 where and Bg are matrices with a dimension of (if x P ) defined by A Q = 

[(cti®ar), (a p ®a p *)] and = [(ai®ai), (a p ®a p )], C s '(a^i^2,T 3 ) is a 

diagonal matrix sized (P x P) defined by C s f (a,xi, 12,13) = diag[ci*(a,xi,x 2> x 3 ),..., 

c/(a,n,T2,X3)] and where c/(ct,xi,X2,X3) is defined by : 
^'(a ,x,,x 2 ,x 3 ) = <Cum(sp(t), s p (t-Ti)\ s p (t-x 2 )\ s p (t-x 3 )) exp(-j27rcu)> 
10 c f> 2 {a ,X|,x 2 ,x 3 ) = <Cum(s p (t), s p (t-Xi)\ s p {t-x 2 ) , s p (t-x z )) exp(-]27ra/)> 

c,/(a ,x,,x 2 ,x 3 ) = <Cum(5p(/), s p (t-x { ), s p (t-x 2 ) , sp(r-x 3 )) exp(-j27ia/)> (22) 

It can be seen that the classic quadricovariance of (6) also verifies 
that Q x = Q x l (0, 0, 0, 0), c p = c p '(0, 0, 0, 0), C s = C s '(0,0, 0, 0). Recalling that 
Q x [i, j, K J\ is the element [N(i -1) + j, N(k - 1) + Z] of the matrix Q x , the 
15 following is deduced: 

& = A Q C t A Q H (23) 
In stating that Q x [i 9 j, fc, Z] is the element [iV(i -1) + I , JV(fc - 1) + j] of 0 t we 
obtain the matrix Q x which is written as follows: 

Q x = B Q C S B Q H (24) 

20 Whitening step 

The first step of the cyclical method orthonormalizes the columns 
of the matrices A Q or Bq contained in the matrices Q x or Q x of the 
expressions (23)(24). The matrices Q x and Q x are Prank hermitian matrices 
and verify the following eigen-element decomposition: 
25 Q x = E x A x E X H et Q x = E x A x E X H (25) 
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Where A v is the diagonal matrix sized (Px F) of the P non-zero values of Q x 
and E x is the matrix sized {ifxPj of the associated eigenvectors. For a full- 
rank matrix B Q > it can be shown that there is an equivalence assuming that 
that the kurtosis values of the sources have a same sign 8 (e = ±1) and 
5 assuming that the eigenvalues of A, also have a same sign 8. In this 
context, the whitening matrix can be built according to T with a dimension 
of (Px N*): 

f - {A,) 1 * E/ (26) 

where (A v ) 1/2 is the inverse of the square root of A r . From the expressions 
10 (24) and (25), it is deduced that: 

ef Q x f H = T B Q (eC s )B Q "f H = \ P (27) 

This expression shows that the matrix T Bq (8C s ) 1/2 with a dimension 
(Px P) is a unitary matrix U . It is then deduced from this that : 

f B Q = U {sC s ) m (28) 

15 It is recalled that the whitening matrix Tof Q x verifies: 

TA Q = U(sC 5 ) m (29) 

Step of cyclical isolation 

From the expressions (28)(29) and (21), it is deduced that: 
20 ^v , (a,r / .r : ,r ? )=7'a l (a,r / .r ? ,r i ) l" = £/(ffC f )" ,/a C I , (a t T, t T 2 ,T 3 ) (eC s )~ m if 

^(a,r/,r:,ri)=ra 2 (a,r/,r:,r j ) f T = C,)' l/2 C, 2 (a,T„T 2f T 3 ) {eC s )~ m U T 
^Vr/.r;,r,)=f ft-V^^r,) f T = U (e C 5 )" l/2 C s 3 (a,T,,x 2 ,T 3 ) (sC s )~ m U J (30) 
When there are Pi sources verifying c*(a,Ti,T2,T 3 )*0 (1< i < Pi), the 
matrix W/(a y r h r 2) r 3 ) is a Pi<P ranking matrix. Thus, the unitary matrices U 
25 and U with a dimension PXP may be decomposed into two sub-matrices 
with a dimension PxPi and Px(P - Pi) such that: 

C/=[£/i U 2 ] and U = [0, U 2 ] (31) 
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where the matrices U\ and U x are sized PxP\ and t/ 2 and U 2 are sized Px(P- 
Pi). The matrices U\ and U x contain the singular vectors associated with the 
non-zero singular values of W x G (a,T h T 2 ,T 3 ). It is deduced from this that: 
W x \a,T h T 2 ,r 3 )=TQ x \a,T h T 2 ,r 3 ) l" = U } (e C 5 )~ m C s W ,,12,13) (eC s ) m U X H 

5 lV x \a y r h r 2t T 3 )=TQ x \a,T h r 2f T 3 ) f T = U^eC,) 1 " C s \a,x u x 2 ,x 3 ) (eC 5 )" m U X T 

WX*.T h T;.T s )=T Qx(<X,T h T 2 , T 3 ) f J =0 x (e C $ C s 3 (a,T lf T J> T 3 ) {S C s )** U y 1 

(32) 

Where the matrix C/(a,n, 12,13) is a diagonal matrix sized PixPi formed by 
non-zero, diagonal elements c*(a,Ti,T 2 ,T 3 ) of the matrix C s *((x,ti,t 2 ,t 3 ). The 
10 matrix C s = C s '(0,0,0,0) sized P\xPi is formed by corresponding elements a 
(1< i < P\). Thus, after a singular value decomposition of W x (ol,ti,t 2 ,t 3 ), it is 
possible to determine the matrices T\ and T x from singular values 
associated with the non-zero singular values and T 2 and T 2 from singular 
vectors associated with the zero singular values such that: 

15 T x = u x n x H , r 2 =c/ 2 /7 2 H , t x = u x n x H ct f 2 =u 2 n 2 H (33) 

where the matrices are 77i, 77 2 , 77, and 77 2 are unitary. From 
W x z \a\n\T 2 \T 3 ) i it is possible to build a matrix Wf{a\ t{ } t 2 ',t 3 ) depending 
solely on the sources of cyclical parameters (a,Ti,T 2 ,T 3 ,e) such that 
c, £ (a,n,T 2 ,T 3 )^0. To do this, the following computation is made: 
20 Wl (a , r, r 2 \ r,')= T? W x \a, r, \ t 2 \ r/) T x = 77, (e C s )~ m C s '(a. r, \ r 2 \ t 3 ) {eC s )~ m 77, H 

W ? ;(a',r/,r ? ;r;)=r l H ^W/' : *,r/) 7^=77, {eC s )~ m C/(a,r/,r/,r/) (<rC 5 )~' /2 77, T 

WfiaW.Tj.Ts)- T x " ^(a,r/,r/,r/) f x *=fi x (s Cy 2 C/(a\T } \T 2 \r 3 )(sC s )' m TI X T 
(34) 

Similarly, from W x e '(a' y Ti',T 2 ',T3) it is possible to build a matrix 
25 W* (a' y Ti',T2',T 3 ) that does not depend on the sources of cyclical parameters 
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(a,Ti,t2,T3,s) such as c,*(a,Ti,T2,T3)=0: Other sources. To do this, the following 
computation is performed: 

^(a,r/,r/,r/) = T 2 " W x \*\t,\t 2 \t s ) T 2 = I7 2 {s C 5 )~ m< C s \a \t,\ r 2 \ r,) (ffC f )" W /7 2 H 
^ v 2 (a , r, z 2 , r 3 )= T 2 H PF v 2 (a , r, r/) f 2 * = 77 2 (* C s ) m C s 2 (a , r/. r,', r/) (*C 5 )" l/2 77 2 T 

5 ^(a,r/,r/,r/)== f 2 H ^ 3 (a",r/J:j/) T 2 *=17 2 (e C,)"" 2 C a \aW.T 3 \T S ){eC t ) ia /7 2 T 
(35) 

where the matrix C/ (a', n', 12,^3] is a diagonal matrix with dimensions (P- 

Pi)x(P-Pi) formed by diagonal elements cf*'[a\Ti' M T2,T3) such that the 

corresponding elements c*(a,rj, T2,T3) are zero elements. The matrix C s sized 

10 (P-Pi)x(P-Pi) is formed by the corresponding elements a (1< i < Pi). 

In particular, in the first version of this variant, it is possible to 

carry out the cyclical isolation in a-0 and e' =1. Writing W[ti',T2,T3) = 

W^\a\Ti , t T2 ) T 3 ) i C s (Ti\T2\T3)=C s t (a',Ti',T2',T 3 ) and C s (ti\t2,T3)= C /(a',rj \x 2 \x 3 \ 

the equations (34) and (35) become: 

15 W x {ti\t 2 \t3)=T 1 h W x {ti,t 2 \t 3 )T x = n x {eCS m C s {ri i T 2t Ts){eC s ) m m H 
(36) 

WATi\T 2 ,T 3 )=T 2 H W x (T lf T 2 ,T 3 )T 2 = /72(^5j" 1/2 ^(ry',r 2 ;r i )(^)" 1/2 /7 2 H 
(37) 

20 Identification step 

The equations (34) and (36) show that it is possible to identify the 
unitary matrices 77! and 77 associated with the sources of cyclic parameters 
(a,Ti,i 2 ,T3,s) in carrying out the joint SVD (singular value decomposition) of 
the matrices Wf (oJ,tj j ) T 2 j,T 3 j) for \<j<K. Thus, to estimate the left unitary 
25 matrix, the joint diagonalizing of the matrices is performed: 

W x £i (a', r, r/ r/) wf (aW, r/ r 3 J f for l<j<K (38) 
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and to estimate the right unitary matrix, the joint diagonalizing of the 
matrices is performed: 

Wf (of, ti J , T 2 J , t 3 ) H Wf (a 1 , rr\ r 2 { r 3 J ) for \<j<K (39) 

To estimate the unitary matrices 77 2 and 77 2 associated with the 
5 sources not associated with the cyclical parameters (a,Ti,t 2 ,T 3 ,e), the joint 
SVD of the matrices Wf (oJ.ti J,t 2 J,t 3 J) for l<j<K is performed in jointly 
diagonalizing the matrices Wf (aJ, r } J, r 2 >, r 3 J) Wf (oJ } t } J,t 2 J,t 3 j) h and then the 
matrices (aJ, n J, t 2 J } t 3 J) h Wf (aJ, r 7 J, r 2 J, t 3 j) . 

Knowing /7i and 77 2 , from the equation (33), the unitary matrices 
10 U] t U 2 and U are deduced to the nearest permutation matrix, in performing: 
t/i= T\ 77, , U 2 = T 2 n 2 et U=[ U\ U 2 ] (40) 

From the equations (9) and (29), it is possible to deduce the 
matrix A Q to the nearest diagonal matrix and permutation matrix, such 
that: 

15 f U =[b } ...b P ]= E x A x m =A Q (£C s ) m ATl (41) 

where 7* is the pseudo-inverse of the matrix T. Each column, bi (1< / < P), of 
the matrix f U is associated with a vector \i q \c q \ l/2 (a^Oq*), 1< q < P, where 
\x q is a complex scalar value such that |n q | = 1. As a consequence, in 
converting each column b\ of the matrix f U into a matrix Bi sized (N x N) 

20 such that B,[f, j] = bi((i-l)N + j) (1 < i,j< N] t it is deduced that: 

B i =^\c q \ n a q a q H for (\<l,q<P) (42) 

In this context, the direction vector Oq of the cf 1 source is associated with 
the eigenvector of Bi associated with the greatest eigenvalue. 

25 Recapitulation of the first version of the cyclical procedure 

In short, the steps of the first version of the cyclical method are 
summarized here below and are applied to L observations x(lT e ) (1< I < L) of 
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the signals received on the sensors ( T e : sampling period). 
Estimation 

Step-1: The estimation of the matrices Q x and (? r from the L observations 
x(lT e ). The estimation of these matrices will depend on the following 
5 assumptions: 

• Stationary and centered case: empirical estimator used in the 
reference [3]. 

• Cyclostationary and centered case: estimator implemented in the 
reference [10]. 

10 • Cyclostationary and non-centered case: estimator implemented in the 
reference [11]. 
Whitening 

Step 2: The eigen-element decomposition of the estimates of the matrices Q 
x and Q x . From these operations of decomposition, the estimation of the 
15 number of sources Pand use of the Pmain eigenvalues such that: 

Qx^Ex A x Ex H et Q x = E x A v £ r H where A x and A r are diagonal matrices 
containing the P eigenvalues with the highest modulus and E x and £ r are 
the matrices containing the associated eigenvectors. 

Step 3: The building of the whitening matrices: T = (A x )~ ia £* H and f = 
20 (A,f I/2 £ r H 

Step 4: The selection of the cyclical parameters (a, n ,t 2 ,13 ,e ) and the 
estimation of the matrix Qx(ol,ti>T2,T3) from the L observations x{lT e ). The 
estimation of this matrix will depend on the following assumptions on the 
signals: 

25 • Stationary and centered case: empirical estimator used in the 
reference [3]. 
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• Cyclostationary and centered case: estimator implemented in the 
reference [10]. 

• Cyclostationary and non-centered case: estimator implemented in the 
reference [11]. 

5 Step-5: The computation of a matrix Wx*(a,T h T2,T 3 ) of (30) from matrices 
Q/(a,rj,r2,r3), Tand T . After singular value decomposition Wx(a, tj,T2,T3), the 
determining of the unitary matrices T\ and T x associated with the non-zero 
singular values and T% and T 2 associated with the zero singular values. 
Selection 

10 Step-6: The selection of the K triplets of delays (ti*,T2*,x 3 *) where 

111*1 + 112*1 + 1X3*1*0. 

Estimation 

Step-7: The estimation of the K matrices Qx(xi fc ,x 2 *,T3*) of (2). As in the step- 
1 this estimation will depend on the assumptions made on the signal such 
15 as: 

• Stationary and centered case: empirical estimator used in the 
reference [3], 

• Cyclostationary and centered case: estimator implemented in the 
reference [10]. 

20 • Cyclostationary and non-centered case: estimator implemented in the 
reference [11]. 
Identification 

Step-8: The computation of the matrices T\ Qx( T i fc > t 2^,T 3 *) Ti H and the 
estimation of the unitary matrix U\ (associated with the cyclical parameters 
25 (a, ti ,t 2 ,13 ,e )) in carrying out the joint diagonalizing of the ^matrices: 
T\ Qjc(xi' c ,t 2 *,T3*) Ti H . 
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Step-9: The computation of the matrices T 2 Qx[n k ,*2 k ,*C3 k ) T 2 H and the 
estimation of the unitary matrix U2 (associated with the other sources) in 
carrying out the joint diagonalizing of the K matrices T2 Q^\ k ^2 k y iz k ) T 2 H . 
Step-10 : The computation of the unitary matrix Uin performing: U =[Ui I/2] 
5 Step- 11: The computation of fu =[b\...bp] and the building of the matrices 
Bi with a dimension (N x N) from the columns bi of T*U. 

Step- 12: The estimation of the signatures a q (1< q < P) of the P sources in 
applying a decomposition into elements in each matrix Bi 

10 

Recapitulation of the second version of the cyclical procedure 

The steps of the second version of the cyclical FORBIUM method are 
summarized here below and are applied to L observations x(lT e ) (l< I < L) of 
the signals received on the sensors ( T e : sampling period). 
15 Estimation 

Step-1: The estimation of the matrices Q x and Q Jrom the L observations 
x(lT e ). The estimation of these matrices will depend on the following 
assumptions: 

o Stationary and centered case: empirical estimator used in the 
20 reference [3]. 

o Cyclostationary and centered case: estimator implemented in the 
reference [10], 

o Cyclostationary and non-centered case: estimator implemented in the 
reference [11]. 

25 Step-2: The eigen-element decomposition of the matrices Q x and Q x . From 
these operations of decomposition, the estimation of the number of sources 
P and the use of the P main eigenvalues such that: Q^E X A X Ex H and Q x = 
E s A T £ V H where A x and A r are diagonal matrices containing the P eigen 
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values with the highest modulus and E x and E x are the matrices containing 
the P associated eigen vectors. 

Step-3: The building of the whitening matrices: T = (Ax)" 2 £* H and f = 
(AJ" ,/2 f r H 

Step-4: The selection of the cyclical parameters (a, n ,t 2 ,13 >e ) and the 
estimation of the matrix Q/(a,rj,r 2j r3) from the L observations x{lT e ). The 
estimation of this matrix will depend on the following assumptions on the 
signals: 

o Stationary and centered case: empirical estimator used in the 
reference [3]. 

o Cyclostationary and centered case: estimator implemented in the 
reference [10]. 

o Cyclostationary and non-centered case: estimator implemented in the 
reference [11]. 

Step-5: The computation of a matrix U^ e (a,rj,r 2 ,r3) of (30) from the matrices 
Q/(cc, t] s T2, t 3 ) , T and T . After a singular value decomposition of 
W x e (a ) T h T2,T 3 ), the determining of the unitary matrices T\ and 7j associated 
with the non-zero significant values and T2 and f 2 associated with the zero 
singular values. 

Step-6: The selection of K sets of parameters (a fc ,n ,c > T2 fc ,T3 /c ,e A: ). 
Step-7: The estimation of the K matrices Q/ k (a t 2 \t 3 ' c ) of (19). As in the 
step-1 this estimation will depend on the assumptions made on the signal 
such as: 

o Stationary and centered case: empirical estimator used in the 
reference [3]. 

o Cyclostationary and centered case: estimator implemented in the 
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reference [10]. 

• Cyclostationary and non-centered case: estimator implemented in the 
reference [11], 

Step-8: The computation of the matrices T\ Qx k (a fc ,xi fc ,i2 ;c ,T3 ;c ) Ti H and the 
estimation of the unitary matrix U\ or (7, (associated with the cyclical 
parameters (a, n ,12 ,13 ,e )) in carrying out the joint SVD of the K matrices: T 
, g/ k (a^T,^,i 2 fc ,T3 fc ) riH. 

Step-9: The computation of the matrices T 2 Q/ k (a fc ,Ti fc ,T 2 ^,i3 fc ) T 2 H and the 
estimation of the unitary matrix U2 or U 2 (associated with the cyclical 
parameters (a, n ,x 2 ,t 3 ,e )) in carrying out the joint SVD of the K matrices: T 
2 Q/ k (a fc ,Ti^T 2 ^,T 3 ^) T 2 H . 

Step- 10: The computation of the unitary matrix Uin performing: U =[Ui U2] 
Step-11: The computation of fu =[b\...b P ] and the building of the matrices 
Bi with a dimension (N x N) from the columns b\ of T*U. 

Step-12: The estimation of the signatures a q (1< q < P) of the P sources in 
applying a decomposition into elements on each matrix Bi 
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